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Abstract. The joint influence of numerical parameters such as the number of particles N, the gravitational softening length 
e and the time-step At is investigated in the context of galaxy simulations. For isolated galaxy models we have performed a 
convergence study and estimated the numerical parameters ranges for which the relaxed models do not deviate significantly 
from its initial configuration. By fixing N, we calculate the range of the mean interparticle separation A(r) along the disc radius. 
Uniformly spaced values of A are used as e in numerical tests of disc heating. We have found that in the simulations with 
N = 1 310 720 particles A varies by a factor of 6, and the corresponding final Toomre's parameters Q change by only about 5 
per cent. By decreasing N, the A and Q ranges broaden. Large e and small N cause an earlier bar formation. In addition, the 
numerical experiments indicate, that for a given set of parameters the disc heating is smaller with the Plummer softening than 
with the spline softening. For galaxy collision models we have studied the influence of the selected numerical parameters on 
\^ ' the formation of tidally triggered bars in galactic discs and their properties, such as their dimensions, shape, amplitude and 

rotational velocity. Numerical simulations indicate that the properties of the formed bars strongly depend upon the selection of 
N and e. Large values of the gravitational softening parameter and a small number of particles results in the rapid formation 
of a well defined, slowly rotating bar. On the other hand, small values of e produce a small, rapidly rotating disc with tightly 
wound spiral arms, and subsequently a weak bar emerges. We have found that by increasing N, the bar properties converge and 
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the effect of the softening parameter diminishes. Finally, in some cases short spiral arms are observed at the ends of the bar that 
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change periodically from trailing to leading and vice-versa - the wiggle. 
Key words, galaxy interaction - numerical simulation - bar formation 



. 1. Introduction force such as, for example, higher powers of the Plummer 

d ' softening ( Athanassoula et al. 2000) or the spline softening 

One of the long-standing issues in A^-body simulations is how ITT V " „ Jlinorl »u JtT"u — lUnmk to. 

& & 3 ( Hernquist & Katz 1989), among others (Dehnen 2001). The 

a specific selection of numerical parameters influence an out- , c c ~ . . , c , , ., ,. , 

r r degree of the force softening is denned by the gravitational 

come. In astrophysical and cosmological simulations there are a ■ . c , f , , ,. 

r 3 b softening parameter s. Small values of s cause high relaxation 

a number of these parameters that may play an important role, . , . • , .. , .. . . , . t t . u 

F ' rates but give a better resolution of the internal structure of the 

and therefore they have to be carefully selected. For instance, . „ , , , c . ■ , , 

3 3 system. On the other hand, a large softening reduces the relax- 

it is well known that collisionless simulations of galaxies are . , . . ., , , , .. . , 

b ation rate, but errors in the force calculation increase because 
performed with a number of particles N, that is far less than c , . .■ c xl . . c 

r r > of the deviation from the Newtonian force, 

the number of stars in a typical galaxy, hence, it induces ar- 

.■ c • , f ... f .... i £ ii t j . The influence of s on the force calculation error was dis- 
tmcial non-uniformities of the gravitational held. In order to I ■ ■ 1 

smooth the field and to avoid the divergence of the gravitational cussed by| Merrit] (| l99fj ), where he found that its optimal choice 
force when the separation between particles becomes small, is related t0 the minimum error that is characterized by the 
, , . .. , , , AA rp, • • so-called bias and variance. He also gives two empirical ex- 
close encounters between particles should be avoided. This is I " 1 1 1 f — 1 

„ ■ , f. ■ f . .. , pressions to estimate s for the Plummer (1911) and Hernquist 

usually done by softening the force between two particles us- j i ' 1 1 * ' * — * 

,-r- , f r ,i xt . ... , (1990) density profiles. This work was later extended to an- 

mg some modified form of the Newtonian gravitational poten- ' 3 v i 1 

, nl ft • • i , • other configurations (Athanassoula et al. 1998, 2000) and to 

tial. For example, the Plummer softening is widely used since & v j m i i [ ' 

the first simulations of Aarseth ( 1963). Other softening algo- other softenin g kernels 0ehnen| | 2OOl|) . Although these crite- 

rithms exist that give a better approximation to the Newtonian ria 8 ive E that minimizes the force errors, they have been de- 

* rived for rigid Monte-Carlo realizations, each of which, when 

Send offprint requests to: R. F. Gabbasov evolved in time, will drift to its individual equilibrium state. 
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Such criteria allow to estimate only the magnitude of s that 
minimizes the initial perturbation due to the modification of 
gravity, but ignores the further evolution. For this reason it is 
important to check whether these criteria apply for dynamical 
evolving systems such as galaxies. 

Two-body relaxation effects and force calculation errors de- 
pend on s and N, and both of them should be simultaneously 
rninimized. It is, h owever, difficult to measure them separately 
jHernquisI 1 993bll . e.g. they man ifest themselves as a joint er- 
ror in the energy of the particles (Hernauist & Barnes 1990). 

Another source of error stems from an inadequate choos- 
ing of the integration time-step At, which also results in poor 
energy conservation. This latter error is typically influenced by 
the maximum acceleration, which is constrained by the soft- 
ening parameter. In addition, there are truncation and roundoff 
errors, but they are small and can be easily controlled. 

The problem of an adequate selection of numerical param- 
eters (N, s, At) for A^-body simulations has been widely dis- 
cussed in the literature, usually separately. In particular, con- 
cerning the effect of the softening parameter in 2D galaxy sim- 
ulations, it was established that in addition to reduction of re- 
laxation, it plays a stabilizing role, and the meaning of a given 
value of the Q parameter s trongly depends on the value of s 
llMillerll97lt Romeo 1994). Indeed, relaxation can heavily af- 
fect th e results of numerical studies of disc stability (e.g,|\Vhite 
Criteria for choosing s for such systems, based on dy- 
namical requirements and ty pe of softening, were discussed by 
iRomeol dl994i 1 19971 [l998). A 3 D study to estimate the nu- 
merical parameters was made by iHernauistl dl987l Il993bl) in 
which the influence of the softening parameter and the number 
of particles on the force errors were investigated. Also, there 
is a rich discussion concerning the choosing of numerical pa- 
rameters and its effects on cos mological sim ulations and small 
scale structure formation (e.g.. |Power et al.l2 003). On the other 
hand, analytical estimations give severe limits on the reliabil- 
ity of A^-body calculations. For instance, it was shown that the 
orbits of particles d iverge ex ponentially with t ime from their 
original trajectories dGoodman et alJl993tlHavesl2003l) . 

For galaxy simulations, the above mentioned effects influ- 
ence any dynamic process, such as the formation of bars in 
spirals that we study in the present work. This problem in- 
volves physical as well as numerical aspects. Since the 70s, 
several numerical simulations have analysed this problem and 
found that cold, rotationally supported stellar discs are glob- 
ally unstable to bisymmetric distortion s and quick ly evolve into 
well defined spontaneous bars (e.g., Hohl 1971). The forma- 
tion of a bar in unstable disc galaxies was extensively studied 
in both two- and three-dimensional simulations (e.g.. | S ^lIwood| 
198lt ISellwood & Carlberd H984T: ISellwood & Athanassoulal 



200C| lAfhanassoula & Misiriotisl 120021) have shown that the 
bar semi-major axis is more than twice as long as the ob- 
served length , which is close to the exponen tial length-scale 
of the disc fElmegreen & Elmegreenj Il985l) . On the other 
hand, it was found that bar rotation velocities were slowed 
down to less than twice the observed velocities in just a few 
Gyrs JPebattista & Sellwood 2000l lAfhanassoula & Misiriotis 



1986; Athanassoula & Sell woodl 1986I) . It was established 
that the bar formation can be su ppressed by i ntroducing 
high enough velocity d ispersions llToomrel 1 19 64). or by a 
massive spherical halo dOstriker&Pee hlesl ll973h . The 2D- 
simulations have shown that the bar's length is not only 
defined by physical parameters but also strongly depends 
on the magnitude of the softe ning parameter JSellw ood 
Il98lt ISellwood &Carlberglll984. Furthermore, recent 3D- 
simulations of isolated galaxies (Debattista & Sellwoodl ll998. 



2002|. This result coincides with predictions of IWein berg 
(1985) that bars are slowed down by dynamical friction with 
a massive halo. However, it was argued that this slowdown is 
associated with the low halo spatial resolution that diminishes 
the ro tational velocity of the formed bar ( Valen zuela & K lvpin 
200 31). A high resolution s imulation with 20 million particles 
(O' Neill & Du binski 2003) still indicates the bar slowdown, al- 
though it is not so drastic. Moreover, it was found that even in 
a stable galact ic disc a bar may emerge du e to discretene ss of 
the ha lo (Walker, Mihos & Hernauist 1996), and Athanas soulal 
(2002, 2003) showed that these effects are related to the disc's 
angular momentum redistribution and resonances of halo par- 
ticles. 

Numerical experiments performed by other authors 
have shown that bars may also form i n disc galaxies 
which interact with a companion galaxy dNogushil Il987t 
iGerin. Combes & Afhanassoulalll990t ISalolll990l). A svste m- 
atic 2D study of tidally triggered bar formation (Sal oi l 199ol) in 
spiral galaxies has shown that bar formation depends simulta- 
neously on the shape of the rotational curve, the disc-to-halo 
mass ratio, the halo concentration, the strength and the ge- 
ometry of the perturbation. These results have been partially 
confirmed by self-consistent 3D simulations dBarnesI Il992t 
iMiwa & Nogushill998UBerentzen et al]l2004h . but the formed 
bars are found to be slower than in isolated models, indicating 
a possible explana tion for the dich otomy of bar characteristics 
found by Elmegreen & Elmegreen ( 1985). 

Motivated by the above arguments, in this work we investi- 
gate for an isolated galaxy in equilibrium and for the collision 
of two galaxies, the range of the numerical parameters (N, e, 
Af) that permit us to minimize numerical artifacts in simula- 
tions. In addition, we compare the spline gravitational soften- 
ing with the Plummer one. For the collision of two disc galax- 
ies we study the influence of the numerical parameters on the 
formation and characteristics of tidal bars. 

This work is organized as follows: In Sect. |5] we describe 
the initial conditions of the galaxy model and the way in which 
the numerical parameters are computed. We obtain a particular 
range for (N, e, Af), which in Sect. |3]is tested to observe de- 
viations from an initial configuration. Based on a convergence 
study we choose the most suitable values for the numerical pa- 
rameters. Then, in Sect.0]we analyse the influence of the pa- 
rameters on the bar formation after the collision of two spirals. 
Finally, in Sect. |5]we discuss the results and in Sect|6]outline 
our conclusions. 

2. Numerical modelling 

A self-consistent galaxy model is usually constructed with par- 
ticles that are governed by the collisionless Boltzmann equa- 
tion coupled to the Poisson equation for the Newtonian grav- 
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itational potential. As mentioned above, the gravitational po- 
tential is softened and, ideally, one should construct a self- 
consistent model for the potential. For simplicity, a standard 
practice is to use an algorithm based on the exact Newtonian 
potential, but this requires to relax the initial galaxy for sev- 
eral crossing times before it can be used for actua l simulations. 
Using this prescription, we follow iBarnesI dl998l) for the con- 
struction of the initial conditions. In this model the bulge and 
the halo are both non-rotating spherically symmetric systems 
with isotropic velocity dispersions. 

The bulge density profile is given by jHernauisl l990): 



Pb(r) 



M h a b 



1 



2n r(r + «b) 



(1) 



and fo r the halo we use Dehnen's y — density profile dPehnenl 
Il99/Sl) : 



Ph(f) 



3M h a h 
4n (r + fl h ) 4 



(2) 



Despite the high concentration of this halo model, it was widely 
used by Barnes in his galaxy collisio n simula tions. The disc 
profile is assumed exponential (Freema nll97(]l) : 



Pd(r, z) = 



M d a 2 

47TZ0 



r sech z - 
Uo 



(3) 



In equations (II 1 3 > Mb, flb and Mh, flh are the mass and scale- 
length of the bulge and the halo, respectively, and Md, a and 
Zq are the mass, the scale-length and the scale-height of the 
disc, respectively. 

Using the Monte-Carlo technique, the density profiles are 
represented by a system of N equal mass bodies, N be- 
ing the total number of particles. For reducing two-body 
relaxation effects and to avoid the mass segregation (e.g., 
Faro uki & Salpeterlll994h we take equal mass particles. The 
number of particles in each component are assigned in propor- 
tion to their masses: Mb : Md : Mh = 1:3: 16. For each 
component, the mass distribution is truncated at a radius that 
contains 95 per cent of its total mass, that would correspond to 
an infinite radius. 

The disc vertical dispersion cr 2 (r) is found from an isother- 
mal sheet equilibrium condition, and the radial velocity dis- 
persion is cr 2 (r) = 4cr 2 (r). The azimuthal velocity disper- 
sion is calculated using the epicyclic approximation cr^(r) = 
cr 2 {r)K 2 /4Q 2 , where k is the epicyclic frequency and Q. the an- 
gular velocity. The net rotational velocit y of the disc is com- 
puted from the asymmetric drift equation ( Bin nev & Trem aine 
Il987l) . The bulge and halo dispersion velocities are calculated 
using equation fl!4i . see below. This equation was integrated 
numerically, including the masses of all components. Then, 
isotropic velocities are drawn from the Gaussian distrib ution. 
For a detai led descrip t ion of the initial conditions see Barnes 
( 1998) and Hernquist ( 1993a). The parameters for our galaxy 
model are summarized in Tabled where the columns from left 
to right give the mass, the number of particles, the cut-off ra- 
dius and the length scale of each component. We use Barnes's 
model parameters and system of units. The disc's scale-height 
is zo = 0.007. The mass, length and time scales are set to 



Table 1. Parameters of the galaxy model. The units of mass and 
length are 2.2 • 10 11 Mq and 40 kpc, respectively. 



Component 


Mass 


Number of 


Cutoff 


Scale- 






particles 


radius 


length 


Bulge 


0.0625 


0.05JV 


1.5 


0.04168 


Disc 


0.1875 


0.15^ 


0.4 


0.0833 


Halo 


1.0 


0.8^ 


6.0 


0.1 




Fig. 1. Rotation curves of the model. The bulge and halo rota- 
tion curves are shown up to the radius of the disc. 



2.2-10 11 M© = 1,40 kpc = 1 and250Myrs = 1 . In these units, 
the gravitational constant is G — 1 . The half mass radius of the 
galaxy is located at R 1/2 ~ 11 kpc. The total rotational curve 
and rotation curves of each component are shown in Fig. [2 

2.1. Methods 

For the time evolution we use a hierarchical General 
Body Smoothed Particle Hydrodynamics (GBSPH) 
treecode written by one of us (M.A.R.M., see details in 
www.astr o.inin.mx/mar/nagbod y) which is similar to Barnes' 
treecode (Barnes & Hut 1986). We use the Plummer softened 
point-mass potential, and the forces are computed up to 
quadruple terms with a tolerance parameter 9 = 0.75. The 
equations of motion are integrated using a second order 
leap-frog algorithm with a fixed time-step. The SPH part was 
switched off. 

To verify our results, we also performed several runs using 
the serial and paralle l versions of the public GADGET code 
(Springel et al. 2001). The main difference between the codes 
is that for the Newtonian potential, GADGET uses a spline ap- 
proximation. The latter calculates for the given Plummer soft- 
ening the equivalent spline softening h = 2.8e. In order to 
maintain numerical similarity, we have fixed the time-step to 
a single value and set the rest of the parameters equal to those 
used by the GBSPH code. Also, the tree structure was updated 
at each time-step. 
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2.2. Selection of the numerical parameters 

For numerical studies one pursues to have a collisionless, stable 
equilibrium system. But, due to discreteness effects and force 
errors this condition is not perfectly fulfilled. If the numeri- 
cal parameters are accordingly chosen, a system that closely 
matches an initial stable equilibrium during a specific time can 
be achieved. Thus, the selection of the parameters is of primary 
importance in any A^-body simulation. 

In spite of the many works on the subject there is still no 
suitable criterion to define the appropriate magnitude of s for 
a given number of particles. Usually, s is chosen in an ad hoc 
way, or by performing various runs with different s and se- 
lecting the one that gives a better conservation of the total en- 
ergy and angular momentum and/or an initial density profile. 
But this procedure demands a large amount of computational 
time and is not always performed. On the other hand, a com- 
mon practice to verify the results is to rerun the simulation with 
higher N, leaving s and At unchanged. 

In fact, given the density profile and the number of particles 
one can roughly estimate the required values of the gravita- 
tional softening parameter and the time-step. We now proceed 
to analyse the criteria for choosing the numerical parameters 
in order to constrain their possible values for the above galaxy 
model. This subsection has a preparatory purpose for the nu- 
merical simulations of the next section, that should further con- 
fine the parameters. The set of parameters that minimize the 
disc heating and total energy errors will be called "optimal". 

2.2.1 . The number of particles 

It is well known that in collisionless A^-body simulations, the 
number of particles should be as large as possible. With the new 
parallel codes it is now possible to simulate galactic systems 
with AT > 10 6 -10 7 . 

The minimum number of particles N required to sam- 
ple a mass M inside an homogeneous sphere of ra- 
dius R can be estima ted from the relaxation time for 
the softened potential ( Huana. Dubinski & Carlber i Il993t 
lAthanassoula. Vozikis & Lam bert 2001 ): 



^relax _ N 

fcross ~ 8 In (R/s) ' 



(4) 



where f cross (Gp) -1 ^ 2 is the crossing time, and p is the mean 
density of the system. Once s is defined and the f le i ax is de- 
manded to be comparable to the age of the universe, the value 
of N can be found. For example, by requiring f re i ax > 10 10 
years, the halo model given by equation @ with f cross « 1 .5 • 10 9 
years and s = 0.01 will need a minimum N m in ~ 3300, 
which is quite a sma ll number. However, numerical simulations 
lAth anassoula et aDl200ll) show that this condition underesti- 
mates the value of the N for centrally concentrated systems by 
roughly an order of magnitude. 

On the other hand, for realistic galactic disc simulations, 
the vertical structure should be resolved, although the general 
structu re in the plane of the disc can be simulated even when 
s > zo dZotov & Morozo\ll987l) . In order to resolve the verti- 
cal structure of the disc component given by equation 0, one 



may require A(R\/2) < Zo, where A is the mean interparticle sep- 
aration estimated within th e disc's half-mas s radius Ri/z, and zo 
is the vertical scale height ( Hernaui stll987l) . For our model this 
gives N^n ~ 10 4 disc particles. In principle, this number can 
be considered as an acceptable minimum for our galaxy model. 
However, in order to make a more realistic estimation we have 
to further analyse the dependence of N with s and Af . 

2.2.2. The softening parameter 

For a given minimum relaxation time and N, the magnitude of 
s depends on the mass distr ibution of the gravitating system. 
In a series of papers Romeoll l 19981 and references therein) has 
analysed the dispersion relation for a 2D disc and discussed the 
criteria for physically consistent values of the softening param- 
eter. It was found that stability and relaxation impose different 
requirements on s. Although this fact also have implications 
for 3D discs, the question in this case is much more compli- 
cated because the particles' vertical motion has to be taken into 
account. 

For an homogeneous 3D mass distribution a natural choice 
for s is to be equal to the mean interparticle separation A. But 
astrophysical systems are centrally condensed and the mean in- 
terparticle separation is determined by local averages, hence s 
is a local parameter. In practice, however, some authors prefer 
to compute the minimum consta nt e that optimize some force 
accuracy level jHernauisl l987). An inconvenience arises for 
very inhomogeneous systems since the election is not obvious. 

There are several popular criteria that are widely used: 



1. IZotov & Morozovl l l 19871) have suggested that the softening 
parameter of a gravitating disc can be taken proportional to 
the radius of interaction, r mU of the particles in its plane, 
which is given by 



M(r) 



(5) 



where m is the mass of an individual particle and M(r) is 
the total mass of the disc inside the radius r. The field is 
considered smooth if the distance at which the interaction 
force between a test particle and its nearest neighbour is 
equal to the force between the same particle and the rest of 
the particles in the disc. The system is said to be collision- 
less if the softeni ng parameter is chosen such that e > r; nt . 
2. lHernquisllfl98 7) have studied force calculation errors, pro- 
duced by small number of particles and the Plummer soft- 
ening of the gravitational potential. He showed that errors 
decrease weakly with increasing N, and are minimized if 
s < A, where A is the mean interparticle separation evalu- 
ated within the half-mass radius R\/2- For s > A the errors 
increase rapidly. The latter is explained, particularly for 
treecodes, by the fact that the expansion of cluster poten- 
tials is accurate only for pointlike particles, and it fails for 
highly softened particles (Hernquist 1987). Then, Newton's 
third law is violated and linear momentum is badly con- 
served. Thus, this criterion can be used to impose an upper 
limit for s. 
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For a spherical gravitating system the value of s that mini- 
mize s force errors can be estimated from an empirical rela- 
tion Jlvlerril 19961) : 



: AN", 



(6) 



where A and a are parameters that depend on particular 
models. For example, for a Dehnen's y — density profile 
of total mass M = 1.0, scale-length a = 0.1 and the cutoff 
radius R = 2 99.8, their values are A = 0.12 and a = -0.27, 
respectively dAfhanassoula et all20 00). However, as it was 
mentioned in the introduction, this criterion only reduces 
the initial shot noise, but does not guarantee that the soft- 
ening chosen in such way will be valid for a long term evo- 
lution. 

Thus, combining these criteria, it is possible to restrict the 
values of the softening for a gravitating system. Let us consider 
a spherical gravitating system of mass M in steady state equi- 
librium. It can be divided in spherical shells of thickness Ar 
whose volumes are 



A V(r, r + Ar) 



(7) 
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Fig. 2. Mean interparticle separation as function of radius for 
Dehnen's y — density profile (upper curve), and the interac- 
tion radius (lower, dashed curve). The curves are plotted only 
up to R = 18 kpc in order to show the behavior near the centre. 
The arrows mark (a) the value of e a calculated from equation 
(|6j, and (b) A evaluated within the half-mass radius. 



Given the density profile p(r), the mass of the shell is 



A M(r, r + Ar) = p(r) d r 



If the system is represented by a Monte Carlo distribution of 
N particles with equal masses m = M/N, then the number of 
particles in that shell is roughly 



A N(r, r + Ar) 



AM(r, r + Ar) 
M 



N. 



In order to estimate the inferior limit for e, the follow- 
ing procedure can be used. The Monte-Carlo particles are in- 

(8) deed agglomerations of ~ 10 4 — 10 6 M when compared to real 
galaxies. Assuming the particles to be Plummer's spheres with 
masses m and half-mass radii rh, they can be treated as point 
masses if the local mean separation A > Ir^. In order to es- 
timate the value of e, we can use the same opening criterion 
9 that is used in treecodes. The structure of the mass distribu- 
tion is unresolved when Ir^/A < 9. The softening radius of a 

(9) Plummer's sphere is related to rj, as 1 3s » rj,, and thus 



The mean interparticle separation in the shell is then given by e — 9 A/2.6 . 



(11) 



A(r, r + Ar) : 



AV(r, r + Ar) 



AN(r, r + Ar) 



1/3 



(10) 



Thus, for each shell a corresponding softening e(r) <x A(r) 
can be found. However, one cannot use a radially varying grav- 
itational softening in the simula tions, since it w ill lead to en- 
ergy non-conservation JZotov & Morozovll987l) . To avoid this, 
one needs a single value of s to be used as the gravitational 
softening parameter for the whole system. When a value of 
s is estimated at some radius R les , it defines the radius of the 
system within which the structure is poorly resolved and the 
obtained results for this region have to be treated with cau- 
tion. Fig. |3 shows A(r) for the above Dehnen's sphere repre- 
sented by N = 262 144 particles together with the values of (a) 
e a w 0.0041 and (b) A(< R l/2 ) ~ 0.0122, that are marked by ar- 
rows. Alternatively, one may use the interaction radius defined 
by equation Q as an estimate for s. The interaction radius ri nt , 
calculated assuming circular orbits, increases linearly with ra- 
dius from zero to roughly 3.31 at the cutoff radius. Using these 
criteria the possible values of softening are confined inside the 
striped area. 



If A is k nown, a rough va lue of s can be calculated. As 
shown by iHernauisl dl987l) . for treecodes the value of 9 e 
(0.7 - 1) represents a compromise between accuracy and per- 
formance. However, there is an inconsistency in implementing 
the Plummer softening due to force calculation as it would be 
between a point mass a nd a Plummer sp here, instead of be- 
tween two point masses dDver&ldll993l) . With the above as- 
sumption, the latter is weakened, but not completely removed 
because of the infinite extension of the Plummer distribution. 
While this inconsistency should be taken into account in sim- 
ulations of systems such as globular clusters, so me authors ar- 
gue it is irrele vant for galaxy simulations (e.g. Romeo 1998; 
IPehnenlEoOll) . On the other hand, the study of Maves! (120031) 
shows that the adequate softening is in the range s e [0.2 - 1] 
times the mean interparticle separation, in agreement with 
equation il Q and Fig. 13 

For each component (bulge, disc and halo) of the galaxy 
model, the mean interparticle separation can be easily eval- 
uated as a function of radius using equations J7I10> . and s 
from equation (1111 . Once s is defined for each component, it 
can be used as an input softening parameter for codes such 
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Fig. 3. Mean interparticle separation as a function of radius for 
the model galaxy for different N. 



as GADGET that handle individual s for each component. 
Individual gravitational softenings are preferable in simula- 
tions, but the total energ y is less conserve d and Newton's third 
law is strongly violated (Dver & Ip Il993l) . To avoid this prob- 
lem individual softenings should be somehow symmetrized in 
order that particles have the same maximum force. Besides, the 
e found in this way does not include the mass distribution of the 
other components. 

If a single value of s is used for a whole galaxy model, then 
the question becomes more complicated, and some averaged 
over component value have to be chosen. Clearly, the single 
value of e will resolve well only some intermediate region of 
the galaxy model, but at the same time, there will be regions 
subject to numerica l heating and poor re solution. 

As is shown bv lAthanassoula et alJ (|2P00) in a compound 
system the force error basically stems from the densest compo- 
nent and depends on its mass fraction. In a galaxy model, high 
density regions are located at the centre and in the plane of 
the disc, while low density regions are found at the periphery, 
mostly in the external halo. Since we are interested in processes 
occurring in the disc component, this should be well resolved. 
Thus, we compute the particle density over the cylindrical vol- 
ume of the disc, which encompasses all three components of 
the galaxy: 



the shape of these curves will depend on the mass distribution 
of a particular galaxy model. 

2.2.3. The time-step 

The accuracy of an integrator scheme is determined by the 
time-step, At, which is in turn constrained by some cri- 
teria, e.g. the Courant condition. If two particles become 
close to each other, their acceleration increases. In order to 
handle their orbits adequately, the time-step should be re- 
duced. Most iV-body codes implement variable time-steps to 
achieve both a better efficiency and accuracy. These codes 
use for each particle an individual time-step, with criteria 
based on its local velocity and/or a cceleration, see for example 
Snri ngel. Yoshida fc White! J200ll) . However, for codes that use 
a single time-step value for the evolution of the whole system, 
the magnitude of At should be carefully estimated. Because the 
maximum acceleration is constrained by the softening param- 
eter and this is related to the number of particles, in choosing 
the time-step one has to consider the role of both parameters. 

In order to find the value of the time-step which is required 
for the correct integration of the orbits of particles in the shell, 
one may use the Courant condition which requires knowing 
their mean velocity of motion. In a spherical system in equi- 
librium steady state the velocities of the particles are usually 
derived from the Gaussian velocity distribution. The velocity 
dispersions, that are functions of the radius, can be obtained by 
taking th e second m omentu m of th e collisionless Boltzmann 
equation (Bin nev & Tre maine 1987J). For an isotropic velocity 
distribution of a spherically symmetric system, <j r — cr^ — erg, 
with 



1 f , „ GM«r') A , 

= ~T\ P {r ) a — dr > 

Pin J r' z 



(14) 



where p(r) is the density and M(< r') is the mass inside a sphere 
of radius r' . The mean square velocity of particles in a shell of 
thickness Ar, is 



< v 2 (r, r + Ar) >- 3 <cr 2 r (r,r + Ar) >, 



with 



<cr r {r, r + Ar)>- 



1 



AM(r, r + Ar) 



J c^(r)p(r)d 3 r. 



(15) 



(16) 



A N(r, r + Ar) = ^ ANi(r, r + Ar), 



(12) The time required to travel the distance A in that shell is 



with 



A Ni(r, r + Ar) = 



AMj{r, r + Ar) 
M~ 



(13) 



The total A{r) for our galaxy model for different N is shown 
in Fig. |3] The volume of the disc was divided in 50 equal vol- 
ume concentric annuli and the mean separation was calculated 
in each annulus using equations (II Oi and (1121 . It is interesting 
that with increasing N, the variation of A shallows. Of course, 



< At(r, r + Ar)>- 



A(r, r + Ar) 
< v(r, r + Ar) > 



(17) 



The space resolution is defined by the radius R Tes , and the mini- 
mum time-step may be taken as the one estimated at R res . Also, 
one can average the value of < At(r, r + Ar) > over the shells, 
obtaining the corresponding to the mean velocity of the system. 
This average time can be used as an estimate of the time-step 
for the orbit integrator. As a consequence, the central region 
having too short dynamical times will not be resolved. 
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For our galaxy model the maximum velocity is determined 
by the rotation curve (~ 250 kms -1 ) and the velocity aver- 
aged over components varies with radius between 160 and 120 
kms -1 . 

Though in galaxy simulations one uses Newtonian dynam- 
ics, the finite propagation speed of the gravitational interaction 
should be included for a realistic self-consist model. If gravity 
propagates with the speed of light c, then information about a 
particle on the opposite side of the system of radius 7? = 150 
kpc, will be received after At = 2R/c * 10 6 years, or ~ 1 /250 
in our units. This implies a restriction on the minimum allowed 
time-step within which the information is received by all par- 
ticles, and/or the maximum radius of the system that can be 
considered in order to maintain the consistency of the simula- 
tion. However, it can be shown that on galactic scales the error 
introduced by particles exterior to cAt is not greater than errors 
caused by the force approximation. The Newtonian approach 
can either be used when a steady state has been reached, be- 
cause then the number of particles in a volume element main- 
tains roughly the same, or when the system is not too extended 
and contribution from distant particles can be neglected. For a 
non-equilibrium rapidly evolving system, such as a cluster of 
galaxies, retarded gravitational potentials should be used. 

In addition, the range of Af is restricted from a numerical 
point of view, by the requirement of stability of the integrator 
algorithm and the required accuracy. 

Thus, the expressions given by equations JlOi and dl7l per- 
mit us to define suitable values for the mean interparticle sep- 
aration and the time-step as a function of the system's radius. 
The advantage is that they can be calculated during the con- 
struction of the model with the number of particles as an input 
parameter. It is clear that the selection of a single value for 
e and At is not obvious. In what follows we will investigate 
the relation between the triad of parameters (N, s, At) and how 
their particular values affect the relaxation of the galaxy model. 



3. Numerical tests of galaxy relaxation 

In order to investigate the relaxation and stability against bar 
formation for various values of s, At, and N, the isolated mod- 
els were evolved up to t = 12.0, which corresponds to 3 Gyrs. 
The main results are summarized in Table fj] where the first 
column gives the model name, the second the total number of 
particles, the third the gravitational softening length, and the 
fourth contains the time-step. The rest of the columns represent 
our results, in the form of control parameters that are described 
below. There are eight series of runs. Models A01-A06 are per- 
formed with different s values that were chosen from Fig.|5]at 
uniform intervals. Models A07-A10 and A16-A23 are similar 
to the first series but for different N. The effect of the time-step 
is studied in models A11-A15. The time-steps are decreased 
by a factor of two, covering the range of velocities 0.64-10 
that satisfy the Courant condition for s = 0.01. Finally, models 
A24- A3 1 are performed using GADGET with different spline 
softenings in order to compare two forms for the force calcula- 
tion. 



The main parameter we use to measure the heating rate of 
the disc is the relative change of the velocity dispersion compo- 
nents. The averaged change in the velocity dispersion is defined 
over n annulus as: 



Yi 



1 _ 

-7 



(18) 



th component of the veloc- 
annulus at time ti and t\, 



where <r\ t (j) and 0^(0 are the k [ 
ity dispersions averaged over the / th 
respectively, and k = [r, <p, z]. Their values are presented in 
columns 5-7 of Table |2 In absence of the heating process, 
these quantities should all be near zero, while non-zero val- 
ues indicate an increment in particle velocity dispersions. The 
system should be relaxed towards a new equilibrium, because 
the initial models are only in approximate equilibrium and be- 
cause of the introduction of the softening. For this we exclude 
the transient time interval t < 0.5 Gyrs. We have evaluated cr 2 
at times t\ = 0.5 Gyrs and tz = 3 Gyrs which gives the ac- 
cumulative estimation of the disc's anisotropic "temperature". 
Accordingly, the values of the relative change of the velocity 
dispersion were averaged over n — 16 annuli. The value of the 
softening that minimizes these quantities will be taken as opti- 
mal. 

For all test runs some disc angular momentum loss is ob- 
served. The transfer rate of angular momentum from the disc 
to the spherical components is related with the efficiency of dy- 
namical friction. For a collisionless system in equilibrium, re- 
distribution of angular momentum is undesired and should be 
minimized. The eighth column of the table shows, in percent- 
age, the relative change of angular momentum with respect to 
the initial value, Z^o- 

The stability of an infinitely thin rotating stellar disc to lo- 
cal axisymmetric perturbations is characterized by the Toomre 
paramete r Q = cr r A c/3.36G£, where S is the disc surface den- 
sity jToomrel fl964). Discs with Q < 1 are found to be sub- 
ject to spontaneous bar formation, otherwise axisymmetric per- 
turbations are suppressed. Moreover, non-axisymmetric per- 
turbations in a differentially rotating disc of finite thickness 
can be suppressed if the critical radial velocity dis persion is 
rough ly twice the Toomre critical velocity dispersion dMorozovl 
1981). In addition, the susceptibility of a thin stellar disc to 
swing amplification of the m = 2 mode is chara cterized by 
the X m = rK 2 j2mnG"L parameter jToomrelll98"l . The nec- 
essary condition to prevent the model from sp ontaneous bar 
growth is Q > 2 JXthanasso ula & Sellwoodl 19861) . and to pre- 
vent swing a mplification in galaxies with steep rotation curves, 
iMihos etalJdl997l) found that X2 > 3 is required. Their values 
in our initial models have a minimum of 1.65 and 3.25, respec- 
tively, and their final values for each run are shown in columns 
9 th andlO th ofTable|2] 

In order to detect the presence of a bar we use the so-c alled 
distor tion parameter, defined as ( Shibat a. Karino & Erig ushi 
2003): 



(19) 
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Table 2. Numerical parameters of galaxy test runs. Tabulated are the model label (1), and the input parameters: the number of 
particles (2), the softening parameter s (3), and the time step At (4). The units of length and time are 40 kpc and 250 Myrs, 
respectively. As a result of simulations the following control parameters are displayed: the relative change of components of the 
disc velocity dispersions (5-7), the disc angular momentum loss (8), the Toomre's Q parameter (9), the Toomre's X parameter 
(10), the average value of the distortion parameter (1 1), the conservation of the total energy (12) and the total angular momentum 
(13) of the system. Finally, we mention the code used (14). 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


Model 


N 


6 


At 




7<p 


7z 


A/ , 

Am 


Q 


X 2 


rj 






Code 


A01 


40960 


0.025 


1/512 


0.782 


0.452 


0.815 


1.5 


2.0 


2.0 


0.065 


0.01 


0.13 


GBSPH 


A02 


- 


0.020 


1/512 


0.670 


0.623 


1.115 


2.1 


2.1 


2.3 


0.033 


0.01 


0.05 


- 


A03 


- 


0.015 


1/512 


0.806 


0.792 


1.089 


2.4 


2.5 


1.5 


0.036 


0.01 


0.12 


- 


A04 


- 


0.010 


1/512 


0.824 


0.605 


1.131 


3.2 


2.5 


1.0 


0.029 


0.03 


0.19 


- 


A05 


- 


0.005 


1/512 


0.863 


0.676 


1.249 


4.5 


2.8 


0.9 


0.026 


0.04 


0.19 


- 


A06 


- 


0.001 


1/512 


0.817 


0.866 


1.374 


6.4 


3.0 


0.8 


0.025 


2.13 


0.12 


- 


A07 


163 840 


0.015 


1/512 


0.457 


0.335 


0.446 


0.8 


1.9 


2.5 


0.021 


0.01 


0.02 


- 


A08 


- 


0.010 


1/512 


0.541 


0.443 


0.532 


1.0 


2.2 


2.5 


0.019 


0.01 


0.02 


- 


A09 


- 


0.005 


1/512 


0.573 


0.417 


0.662 


1.3 


2.3 


2.0 


0.015 


0.02 


0.05 


- 


A10 


- 


0.001 


1/512 


0.679 


0.556 


0.909 


2.8 


2.6 


1.5 


0.016 


0.57 


0.03 


- 


All 


- 


0.010 


1/64 


0.562 


0.398 


0.477 


1.2 


2.1 


2.0 


0.019 


1.08 


0.07 


- 


A12 


- 


0.010 


1/128 


0.527 


0.397 


0.466 


1.0 


2.1 


2.0 


0.017 


0.04 


0.05 


- 


A13 


- 


0.010 


1/256 


0.536 


0.462 


0.470 


1.1 


2.2 


2.5 


0.019 


0.02 


0.02 


- 


A14 


- 


0.010 


1/512 


0.541 


0.443 


0.532 


1.0 


2.2 


2.5 


0.019 


0.01 


0.02 


- 


A15 


- 


0.010 


1/1024 


0.523 


0.482 


0.478 


1.2 


2.2 


2.2 


0.025 


0.02 


0.05 


- 


A16 


491 520 


0.010 


1/512 


0.347 


0.269 


0.195 


0.6 


1.9 


2.6 


0.016 


0.01 


0.04 


- 


A17 


- 


0.008 


1/512 


0.386 


0.224 


0.263 


0.5 


2.0 


2.8 


0.012 


0.01 


0.02 


- 


A18 


- 


0.006 


1/512 


0.384 


0.270 


0.277 


0.6 


2.0 


2.8 


0.011 


0.01 


0.02 


- 


A19 


- 


0.004 


1/512 


0.340 


0.277 


0.316 


0.7 


2.1 


2.8 


0.009 


0.02 


0.02 


- 


A20 




0.002 


1/512 


0.333 


0.253 


0.382 


0.7 


2.2 


2.7 


0.008 


0.04 


0.03 




A21 




0.001 


1/512 


0.315 


0.256 


0.455 


0.8 


2.2 


2.6 


0.008 


0.21 


0.03 




A22 


655 360 


0.008 


1/512 


0.320 


0.247 


0.155 


0.4 


1.9 


2.8 


0.010 


0.02 


0.01 




A23 




0.001 


1/512 


0.306 


0.253 


0.389 


0.7 


2.1 


2.5 


0.009 


0.12 


0.02 




A24 


163 840 


0.015 


1/512 


0.786 


0.618 


0.536 


2.0 


2.4 


2.4 


0.035 


0.09 


0.02 


GADGET 


A25 




0.010 


1/512 


0.861 


0.648 


0.566 


2.6 


2.5 


2.1 


0.030 


0.12 


0.06 




A26 




0.005 


1/512 


0.600 


0.477 


0.708 


2.0 


2.5 


1.9 


0.016 


0.13 


0.05 




All 




0.001 


1/512 


0.658 


0.529 


0.908 


2.5 


2.6 


1.4 


0.014 


0.79 


0.05 




A28 


655 360 


0.008 


1/512 


0.674 


0.505 


0.271 


1.3 


2.1 


2.6 


0.018 


0.21 


0.57 




A29 




0.001 


1/512 


0.308 


0.276 


0.425 


0.6 


2.2 


2.7 


0.008 


0.30 


0.42 




A30 


1310 720 


0.006 


1/512 


0.369 


0.276 


0.141 


0.3 


2.0 


3.1 


0.010 


0.11 


0.29 




A31 




0.001 


1/512 


0.267 


0.229 


0.312 


0.5 


2.1 


2.8 


0.008 


0.13 


1.11 





where 

I - 1 21 

'xx l yy ^xy 

^ = nT' = T^TT' (20) 

*xx ~ 1 yy l xx ~ 1 yy 

and the moments of inertia are given by 

hi = 2^ m k x k x k > *> J - (*> y) ■ ( 2 1 ) 

k=l 

This parameter allows us to detect any non-axisymmetric de- 
formations such as bars in the disc plane. The mean values of 
77 for each run are summarized in the 1 1 th column of the same 
table. We have found traces of a bar for 77 <: 0.02, which repre- 
sents a rough threshold. 

As criteria to estimate global errors in the force calculation, 
we use the conservation of the total energy, E, and the total an- 
gular momentum, L, whose percentage relative errors are also 
displayed in Table |2] (columns 12 and 13, respectively). 



In general, the models are quite stable and become dynam- 
ically "hot" (Q > 2) which should prevent them from sponta- 
neous bar formation. However, at the same time X2 is reduced, 
making the disc susceptible to bar mode amplification. The ma- 
jor disc heating indeed occurs for t < 0.25 Gyrs, while the sys- 
tem shifts to a new equilibrium; later on, the velocity disper- 
sions grow relatively slowly. From Table|3]a general tendency 
of Q to decrease linearly with increasing s can be observed for 
the Plummer softening. 

As can be seen from Table |5] in models A01-A10, the vari- 
ation of s affects most of the control parameters. There is a 
clear tendency of the parameters y r , y^ and y z to increase with 
decreasing s, caused by the heating process of the disc. The 
tendency is pronounced for small N, and weakens for large N 
(models A16-A23). Indeed, no simple interpretation of the be- 
haviour of y r and y$ on e can be drawn for the simulations 
with higher number of particles (models A16-A23) and those 
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performed with GADGET (models A24-A31). However, it is 
clear that with the increasing N the heating is reduced. 

During the evolution, there are some halo particles that es- 
cape from the system, but their number is relatively small (< 2 
per cent), and for all models the density profile of the spherical 
components is well preserved. For small values of e, the disc 
rapidly heats up, and reaches very large values of Q (> 2). In 
these cases, we observe suppression of a transient spiral pattern 
after a few rotation periods, indicating the presence of dynam- 
ical heating of the disc. For large e the central region of the 
galaxy model with Af = 40 960 rapidly collapses, and a ring- 
like structure moving outwards is observed during the first few 
time steps. The latter is due to the non-equilibrium evolution 
of the model caused by the lowered potential energy and the 
strong deviation from the Newtonian law. Regarding the disc 
vertical structure, it is observed that for smaller s, discs become 
thicker during the evolution. The two-body relaxation effects in 
the vertical direction become important as the vertical structure 
is being resolved. This makes the disc hotter and thicker (com- 
pare y z values in Tabled. 

We allowed some models to run up to t — 8 Gyrs, at which 
77 reaches the value ~ 0.1. Since t « 4 Gyrs we found a weak 
diffuse bar, that appears even for relatively high values for Q 
(> 2.0), that indicates that such high Q-values alone cannot 
suppress the bar instability. A more massive hal o also prevents 
the ba r formation for a longer time. However, lAthanassoulal 
has recently shown that a live halo can play a rather 
destabilizing role. The time at which the bar forms varies 
slightly with Af and e, appearing earlier for small Af and large 
s (see values of 77 of the Table |2j. This confirms findings of 
IWalker et alJ dl996l) . who showed that the bar emerges later 
when Af is increased, and for a rigid halo the bar does not form. 
Large softenings lead to formation of a bar that is more pro- 
nounced, and appears at earlier times (t « 2 Gyrs) in spline 
softening models than in Plummer's. 

To demonstrate that the total energy conservation mainly 
depends on the time-step, we vary At and fix the rest of the nu- 
merical parameters constant. This is done in models Al 1-A15, 
which show that the main change among the control parame- 
ters is in the energy. The deterioration of energy conservation 
for smaller gravitational softenings is explained by the errors 
introduced by the numerical leap-frog integrator for particles 
with high accelerations which are integrated with large time- 
steps. When the softening is small but the time-step is suffi- 
ciently large, such as for example in model A06, a systematic 
non-conservation of the total energy is observed at each time- 
step. As shown in models A11-A15, for s = 0.01 and for typ- 
ical velocities v ~ 1 or ~ 156 kms~', a time-step of 1/128 
is quite sufficient to obtain energy conservation of ~ 0.04 per 
cent after 3 Gyrs of evolution. On the other hand, comparing 
the energy conservation between models with e = 0.001, it can 
be noted that the increase of Af also improves energy conserva- 
tion. This is probably due to the accomplishment of a smoother 
gravitational field. 

For comparison, we have performed several runs using the 
GADGET code (models A24-A31). We have decided to skip 
the runs with Af = 491 520 but preformed instead the simula- 
tions with N = 1 310720. These runs, in general, show higher 




1,00 



,0,99 



0,98 



0,97 

0,0 0,5 1,0 1,5 2,0 2,5 3,0 
time, Gyrs 

Fig. 4. Evolution of the disc angular momentum for models 
A24-A27, normalized to the initial value. 



heating rates and poorer energy conservation. The higher heat- 
ing rate of these models suggests that the Plummer softening 
makes the system less collisional than the corresponding spline 
softening. This is due to a slo wer converge nce of the Plummer 
model to the Newtonian force llTheisll998h . 

The high particle numbers of models A16-A21 and A22- 
A23 decrease the overall process of dynamical heating and di- 
minish the distortions of the disc surface density. By comparing 
the y values with the same s and different number of particles, 
we found that augmenting N by four times, the y values dimin- 
ish few times. For the given range of s in models A16-A21, the 
weak variation of most of their control parameters indicates 
that they are close to convergence. Models A22-A23 confirm 
this convergency. Models A30-A3 1 that were performed using 
GADGET code also show further reduction of the disc heat- 
ing and give results comparable to those of models A22-A23, 
which were performed with GBSPH code. Note also, that these 
models show no signs of a bar. 

The effect of the softening clearly manifests itself in the 
amount of angular momentum transferred from the disc to the 
bulge and halo. The eighth column of Table|2]shows the depen- 
dence of the disc's angular momentum loss on the numerical 
parameters. In all runs decreasing the Plummer softening, in- 
creases the disc angular momentum loss. The observed loss rate 
is almost constant through a single simulation and is different 
for each run. For spline softening (models A24-A31) this pro- 
cess is slightly different, as can be seen in Fig.|4] The loss rate is 
linear in models A26 and A27 as for models with the Plummer 
softening, but in models A24 and A25, as well as in A28, the 
bars emerge quite rapidly at t « 2 Gyrs, which increases the 
growth rate. Concerning the total angular momentum, it can be 
noted that it is conserved quite well and almost independent on 
our chosen numerical parameters. 

The numerical tests have shown that for the models in 
Table |2 the softening s = 0.001, for the considered range 
of N and Af, has the higher heating rate and the worst to- 
tal energy conservation, but, at the same time, smaller non- 
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axisymmetric distortions. This value probably fulfills a criteria 
similar to the ones given through equations (0 and for a 
compound galaxy model and minimize the initial force errors. 
However, these models cannot be accounted as realistic ones 
(except for the model A31), but, for comparison, we will use 
them in some of the collisions described in the next section. 
Other extreme cases are those that use e evaluated at the edge 
of the disc. Despite the fact that in some models they do not 
resolve the inner and vertical structure of the disc, but only its 
outer part, they show the smallest heating and a good conserva- 
tion of E and L. These models are also used in collision runs. 
From Table|5]we choose as compromise values for e opt the ones 
that correspond to models A03, A08 and A 19. In addition, the 
models A22-A23 and A28-A29 have also been taken since they 
show convergence in the control parameters. 

4. Galaxy collisions 

In this section we describe the effect of the number of particles 
and the gravitational softening on the properties of an interact- 
ing system. Galaxy models are constructed with different e and 
N values and used in galaxy collision simulations to study the 
effect of s on the formation of bars and their properties, such as 
the amplitude and rotational velocity. In order to reduce the pa- 
rameters' space, all collision simulations were performed with 
the single time-step At = 1 /256. 

For close galaxy encounters, interaction between particles 
are much more important than for an isolated galaxy. The value 
of e opt found for an isolated steady state system is not neces- 
sarily optimal for an interacting system where strong non-linear 
effects and large density gradients are developed. Force errors 
and relaxation processes are drastically increased by the en- 
counter due to the almost head-on collision of particles along 
the shock front. Moreover, the time-step found for the isolated 
galaxy models may be inadequate for collision simulations. 

For the simulations, equal galaxies are placed on parabolic 
orbits, calculated from the two-body problem with time to peri- 
centre f p = 0.75 Gyrs and pericentric separation p = 2Rd = 32 
kpc. The latter is chosen such that at the pericentre the discs are 
just touching each other, and the main perturbation is provided 
by the interaction of the extended haloes. This pericentric sep- 
aration provides enough time for the bars to evolve and, at the 
same time, to investigate the merger process. With such orbital 
parameters, the galaxies are initially separated by R p as 160 
kpc, which sets the system with overlapped haloes and thus, it 
introduces some initial perturbations. The relatively small R p is 
chosen to reduce cumulative numerical errors, where the pre- 
ferred duration of interaction should be as small as possible. 
The galaxies were relaxed up to t — 0.25 Gyrs before placed on 
their orbits. 

We only consider prograde-retrograde collisions that allows 
us to investigate, with a single configuration, bar formation pro- 
cesses in the two possible directions of rotation. Planar col- 
lisions (i.e., disc planes of the galaxies coincide) are chosen, 
because this configuration is the most violent and provides the 
maximum perturbation and maximum transfer of orbital to in- 
ternal angular momentum of the discs. In the present study we 
do not pretend to cover the entire orbital parameters space. The 



Table 3. Parameters of collisions. The units of time and length 
are 250 Myrs and 40 kpc, respectively. 



Model 


N 


At 


E 




i-0 


Code 


BOl 

B02 
B03 


40 960 


1/256 
1/256 
1/256 


0.001 
0.015 
0.025 


14.53 
0.09 
0.14 


0.06 
0.10 
0.40 


GBSPH 
- 


B04 
B05 
B06 


163 840 


1/256 
1/256 
1/256 


0.001 
0.010 
0.015 


5.40 
0.11 
0.10 


0.05 
0.10 
0.16 


- 


B07 
B08 
B09 


491 520 


1/256 
1/256 
1/256 


0.001 
0.004 
0.010 


2.17 
0.10 
0.09 


0.07 
0.16 
0.11 


- 


B10 
Bll 


655 360 


1/256 
1/256 


0.001 
0.008 


1.61 
0.07 


0.06 
0.08 




B12 
B13 
B14 


163 840 


1/256 
1/256 
1/256 


0.001 
0.010 
0.015 


7.0 
0.33 
0.27 


4.98 
5.87 
4.98 


GADGET 


B15 
B16 


655 360 


1/256 
1/256 


0.001 
0.008 


2.17 
0.11 


2.94 
2.33 





emphasis is rather to investigate the influence of the numerical 
parameters on the bar formation for a given configuration. 

Details of collision models and conservation of total en- 
ergy and angular momentum are summarized in Table [3] 
Animations of some collision simulations are available at 
www. astro, inin. mx/ruslan/tidalJbars. 

In order to perform a reliable analysis of the bar charac- 
teristics, its centre of mass should be well defined. The disc 
particles that are ejected to large distances can spuriously af- 
fect the parameters under study, and therefore these particles 
are excluded from the analysis. We then calculate the centre of 
the disc as follows: First we compute the centre of mass using 
all the particles of the disc. We then recalculate a new centre 
of mass with only the disc particles that are located inside 2Rd 
with respect to the firstly calculated centre of mass. The centre 
of mass found in such a way gives roughly the position of the 
centre of mass of the disc with respect to which the bar parame- 
ters are further computed. However, this procedure fails during 
the last stage of the collision process when both systems are 
found inside the same radius. 

To characterize quantitatively the bar, we use two criteria. 
The first is given through the distortion parameter rj defined 
above. The second is to compute the harmonic amplitudes of 
the disc mass distribution in the equatorial plane, which indi- 
cate the presence of any non-axisymmetric deviation. The nor- 
malized coefficients are (Sell wood & At hanasso uiall986h 

1 - 

A m = - expO'mfy), (22) 

where n is the number of particles within the radius of the disc 
Ri, and Gj is the polar angle of the j' h particle with respect to the 
centre of mass. Since we are mostly interested in bisymmetric 
perturbations, we set m = 2, but we also investigate ampli- 
tudes up to m = 8 to compare among the harmonics. We have 
found that the 77 and A2I bar strengths criteria give similar re- 
sults. However, stronger fluctuations are observed using 77, and 
in what follows we will only refer to IA2I. Some morpholog- 
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Fig. 5. The separation between the centres of the two discs as a 
function of time for models B04-B06. 




0,0 

012345678 



time, Gyrs 

Fig. 6. The evolution of the amplitude of the second harmonic 
for models B01-B03. The solid line corresponds to the first 
galaxy (prograde orbit) whereas the dotted line corresponds to 
the second galaxy (retrograde orbit). 



ical characteristics of the bar can be obtained from projected 
density contours and projections of particles on the orthogonal 
planes. 

An important quantity that characterizes the bar is its ro- 
tational velocity Q, that is given by the rate of change of the 




time, Gyrs 



Fig. 7. The evolution of the amplitude of the second harmonic 
for models B04-B06. The correspondence of curves is the same 
as for the previous figure. 



angle between the coordinate axis and the principal axis of the 
moment of inertia of the bar. We divide the plane of the disc 
into 16 concentric cylinders for which the eigenvalue of the 
moment of inertia tensor I xx is calculated. Then we iteratively 
rotate the coordinate system with an angular step of 2° until the 
maximum value of I xx is reached for each annuli. The corre- 
sponding averaged polar angle of rotation with respect to the 
original coordinate system is taken as the phase of the bar. The 
procedure is repeated for each snapshot file, and O is found 
as d(pldt, where dt = 1/32 is the time interval between snap- 
shots. To correctly identify the bar major axis, we impose the 
restriction on the minimum value r] + = 0.15, which excludes 
the time intervals when the bar is not present and excludes also 
the central annulus, which is almost axisymmetric. The second 
restriction consists to discard annuli with a phase difference of 
more that 10° from the averaged phase of the inner annuli. This 
procedure excludes the annuli where spiral arms and the ring 
are located. 

For all galaxy collision models, the first encounter occurs at 
time t « 0.8 Gyrs and the second at t « 6 Gyrs. The tidal inter- 
action induces a bipolar perturbation that develops spiral arms 
in the first galaxy which is in prograde orbit, while the second 
galaxy does not show any significant deformation. The subse- 
quent evolution is different for each model as will be described 
below. 
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Fig. 8. The evolution of the amplitude of the second harmonic Fig. 9. The evolution of bar angular velocity for collision runs 
for models B07-B09. The correspondence of curves is the same B01-B09. Shown only £2 of the first (direct) disc, 
as for Fig. [6] 



Fig. [5] shows the separation between the centres of the two 
discs as a function of time for models B04-B06. It can be noted 
that in model B04, after the first encounter, the discs have a 
maximum separation distance at t « 3.5 Gyrs, which is less 
than those of models B05 and B06 due to a smaller soften- 
ing of the force. In addition, the whole collision process is 
slightly faster for model B04, and slower for B05 and B06. 
By analysing the curves for runs B01-B03 (not shown here) we 
have noted that they have a similar shape, and that by increas- 
ing N (collisions B07-B09) the difference tends to vanish. 

Models B01-B09 were performed with the GBSPH code 
with three different values of N and s. We first describe col- 
lision runs B01-B03 which were performed with N = 40 960 
particles per galaxy. As can be seen from Table[5] the run B01 
shows bad conservation of the total energy due to integration 
errors caused by inadequate time-step for such small s. But at 
the same time, there is a good angular momentum conservation. 
The evolution of the bar amplitude for these runs is shown in 
Fig lUfor both prograde (solid line) and retrograde (dotted line) 
discs. There are also small contributions from other harmonics, 
mostly from even ones. 

After the first encounter, transient spiral arms are formed in 
the prograde galaxy which then become more and more tightly 
wound until they form a disc again. At the same time, a small 
diffuse bar emerges which maintains an oval shape of apparent 
length / ~ 10 kpc until the second collision, after which the bar 



is amplified to length Z — 16 kpc and develops a butterfly shape 
when viewed from the disc's plane. 

In model B03 after the first passage, transient open spiral 
pattern is generated in the first galaxy, and at the same time, 
a strong bar of length I ~ 20 kpc appears. Note that this is 
the simulation with the longest bar. The spiral arms begin to 
wind around the bar, forming a ring around it which is then 
slowly absorbed by the bar. The bar maintains its length until 
the second close approach, after which it is partially destroyed 
and reduced to / ~ 16 kpc. 

The second galaxy that is in retrograde motion shows a 
slowly growing bar mode, which is amplified by the second en- 
counter (see Figs. I6I8> . The barred galaxies collide again and 
finish in a disc-like remnant. The initial growing mode of the 
second galaxy that becomes notable after t « 4 Gyrs is sim- 
ilar to the one in the isolated evolution, described in sec. [3] 
Therefore, it seems not to be caused by the tidal interaction, 
but instead by the instability of the galaxy model. This is ob- 
served in all collision models presented in Table [3] The insta- 
bility growth rate is slow and by the time of the first encounter 
the bar has not yet been formed. However, the evolution of the 
bar amplitude after ( s 6 Gyrs may be influenced by both ef- 
fects. 

The behaviour of the B02 collision run, is similar to the 
B03, except that the bar is shorter (/ ~ 16 kpc). After the second 
collision, the bar is strongly amplified, because at the pericenter 
it becomes aligned with a diffuse bar of the second galaxy. The 
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Fig. 10. Projected density contours of the first (prograde) disc 
in units of 137.5 Mq/pc 2 . Columns from left to right corre- 
spond to models B04, B05 and B06. The boxes are 32 x 32 kpc 
long. Contours are shown for times (panels from top to bottom) 
0, 1.5, 3.5, 5.5, 6.5 and 8 Gyrs. 
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Fig. 11. Projected density contours of the second (retrograde) 
disc. Columns from left to right correspond to models B04, B05 
and B06 respectively. All scales are the same as in the previous 
figure. 



amplitude of the second harmonics of the bar reaches the value 
~ 0.6, but apparent length of the bar does not change, obtaining 
rather a butterfly shape. 

The models B04-B06 are performed with N = 163 840 par- 
ticles per galaxy and show a behaviour similar to models B01- 
B03 but with richer details in the inner structure. For example, 
a §- shaped structure surrounding the bar was present in these 
runs, but it was not observed in simulations with N = 40 960 



particles. The notable tidal perturbation decay is also observed 
in model B04, but the difference between the models B05 and 
B06 is not so strong (Fig. 0. The models B07-B09 are per- 
formed with higher N and show a more detailed inner structure, 
but only a weak dependence of bar strength on s (see Fig. [HJ. 

Similar differences are observed in the rotational velocities 
of the bars. The evolution of bar angular velocities for mod- 
els B01-B09 are plotted in Fig. [9] As it can be noted form the 
upper panel, the magnitude of Q differs significantly among 
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Fig. 12. The evolution of the amplitude of the second harmonic 
for models B12-B 14. The correspondence of curves is the same 
as in the Fig. [6] 



the runs. However, with the increase of N these differences are 
reduced, and the shape of the curves becomes similar (mod- 
els B07-B09). A clear evidence for correlation between the bar 
strength and its angular velocity can be seen when Figs. l6l9l are 
compared. For example, after the second encounter the decay 
of £2, observed in models B07-B09 occurs when the bar, due 
to alignment with the second bar, is amplified. This reduction 
can be associated with the angular momentum exchange which 
amplifies the bar and reduces its patt ern speed, in accordance 
with results of Berentze n et alJ (|2004). 

The softening affects not only the shape of the bars, but also 
its internal structure. Figs. ^|and^2p resen t the contours of 
the projected density for models B04-B06 at different t. From 
the upper panels it can be seen that due to different s, the initial 
discs have different central densities, which may be responsible 
for their subsequent evolution. The rest of the panels show that 
bars with smaller s are rounder, denser and smaller than those 
with larger e. 

Summarizing, the results of models B01-B09 have shown 
that the decrease of s damps perturbations in the prograde disc. 
For the second galaxy which is retrograde, the amplitude A2I 
slowly grows with time, and the growth rate slightly depends on 
the chosen s. Figures l6l8l indicate that the initial perturbation 
has almost the same peak value of A2I <; 0.5, but further decay 
of tidal response depends on s and N. 




Fig. 13. Projection of particles of the first disc for run B06 for 
times: a) 3.125, b) 3.875, c) 4.375 and d) 5.125 Gyrs. For face- 
on and edge-on projections the boxes are 40 x40 kpc and 8 x 40 
kpc long, respectively. The bar always rotates clockwise, but its 
arms wiggle from trailing to leading and vice-versa. 



For comparison, we have performed collision simulations 
using a parallel version of the GADGET code (runs B12- 
B16). As can be seen from Table[5] the total angular momen- 
tum is poorly conserved. After the first encounter, this non- 
conservation shifts the galaxies from the orbital plane, which is 
more notable for small N (runs B12-B14). As a consequence, 
further encounters are not planar and the merger process takes 
longer. Despite this, the first stage of the bar formation can be 
compared with previous runs. Fig. II 21 shows the evolution of 
the amplitude of the second harmonic for collisions B12-B14, 
which is similar to models B04-B06 except that the amplitude 
of the bar is a bit smaller and a stronger bar develops in the sec- 
ond disc. The collision simulations performed using GBSPH 
and GADGET codes for N = 655 360 show no significant dif- 
ferences in IA2I for the same e; these curves are not shown. 

Finally, we will describe some peculiar characteristics of 
collision and merger processes for different softenings. A gen- 
eral interesting feature is observed between the first and the 
second collisions. The two large spiral arms that are formed 
just after the first collision in the prograde galaxy become more 
and more tightly wound, and form a ring-like structure that sur- 
rounds the diffuse bar. This ring is then slowly destroyed ac- 
creting into the bar, and short transient arms are formed at the 
extremes of the bar (Fig. 1131 . These arms change from trail- 
ing to leading and vice-versa with a period T^m ~ 0.675 Gyrs 
(model B06) though the direction of rotation of the bar does not 
change. We refer to this dynamical feature as the wiggle effect. 
As most particles are absorbed by the remaining pronounced 
bar, the wiggle slowly vanishes but does not completely disap- 
pear and is sustained until the second collision. The wiggle is 
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observed in runs where a clear bar appears and in all runs with 
N = 491 520. In the same figure, shown aside, a classical bar 
buckling leading to the box/peanut shape may be seen. 

In our orbital configuration, after the first encounter up to 5 
per cent of the particles of the first disc leave the initial radius, 
and up to 10 per cent after the second. The second disc looses 
only a few particles during the collision. The particle loss rate 
is almost independent on the numerical parameters and in each 
case the merger remnant is composed of roughly 70 and 85 per 
cent of particles of the first and second discs, respectively. In 
runs B01-B06 the remnants have elongated orbits due to par- 
ticles of the prograde disc, whereas particles of the retrograde 
disc form a symmetric constant density core (Figures ^3 EH 
bottom panel). The mergers with large softenings finish in a 
disc-like remnant, whereas for small softenings, the remnant 
has rather an elliptical structure. 

5. Discussion 

The aim of the convergence study was to investigate how the 
numerical parameters influence the relaxation and stability of 
a galaxy model against bar formation. We first discussed how 
the parameters N, s and At can be restricted, and for a particular 
galaxy model obtained valid ranges for e and At. For a galaxy 
model represented by a given number of particles, the local 
mean interparticle separation in the plane of the disc can be 
easily estimated. Although the gravitational softening parame- 
ter can be chosen to be equal to any particular value of A within 
the established range, one still has to make a compromise be- 
tween force accuracy and relaxation. But it seems that there are 
no universal choices, and the selection of a particular e opt value 
for for a whole system depends on dynamical aspects one wants 
to study ( Romeo 1998). In this context we wish to minimize the 
relaxation effects that are responsible for drifting of a system 
from its initial condition. The convergence study performed in 
section[3]showed that relaxation processes are minimized for e 
equal to the maximum A(Rd), which is evaluated at the edge of 
the disc. This s gives low heating and good energy conserva- 
tion, but it neither resolves the vertical nor the radial disc mass 
distributions. On the other hand, the minimum s = A(0) is un- 
acceptable in simulations because of the high heating rate, and 
because it demands a too small time-step to achieve acceptable 
energy conservation. In Table |2 we have presented a series of 
models with gravitational softenings chosen within the estab- 
lished range which is defined by N. Fig.[5]indicates that the A 
range decreases with increasing N. 

We have found that for our galaxy models the characteris- 
tics of the resulting systems tend to converge for Af ^ 491 520. 
Namely, the reduction of the range of A by increasing N also 
reduces the variation of control parameters. The parameters y z , 
Q and the disc angular momentum loss rate, used to measure 
the disc heating, proved to be quite informative. The disc an- 
gular momentum loss rate gives the efficiency of the dynamical 
friction between the disc and the halo, whereas the increase 
of the vertical velocity dispersion, observed in all runs and for 
both softening methods, measures the thickening of the disc. 
However, it seems that those quantities are related. The dynam- 
ical friction is enhanced by increasing the thickness of the disc, 



and this, in turn, increases the angular momentum loss rate. 
Indeed, the change in the slope of angular momentum loss of 
the disc serves as a good indicator to identify the bar's emer- 
gence. As can be seen from Tabled simulations with 491 520 
particles greatly reduce the disc angular momentum loss, which 
is less than 1 per cent within 3 Gyrs, and only weakly depends 
on e. The disc thickening is also reduced to about three times, 
when compared to runs with 40 960 particles. Simulations with 
N - I 3 10 720 show that the heating process is weakly affected 
by s in comparison with the simulations with smaller particle 
numbers for the corresponding ranges of A. From the same ta- 
ble it can be noted that final values of Toomre's parameter Q 
for the Plummer softening show approximately a linear depen- 
dence on e. 

The heating of the disc can be reduced by increasing either 
the softening parameter or the number of particles. However, 
increasing s within the allowed range reduces the heating in a 
slower way than increasing N (compare the values of y in the 
Table |3. This reaffirms the common belief that the increasing 
of the number of particles lets to effectively decrease the nu- 
merical effects and at the same time to increase the resolution. 
To obtain the minimum N that resolves the disc vertical struc- 
ture, it is necessary that s 5= zo- The criterion of the interaction 
radius, given through equation (0, showed to be inadequate 
because it gives too small s values. On the other hand, the cri- 
terion given through equation dl 1 i would give a minimum e, 
which is approximately 0.3 times the mean interparticle sepa- 
ration, estimated within the half-mass radius. The interaction 
force in this case will be well represented and the bar will ap- 
pear later, but the system will be heated. 

On the other hand, the time-step, which is proportional to 
the gravitational softening parameter, defines the accuracy of 
the integration of the orbits and is the main factor responsi- 
ble for energy conservation. The criterion for the time step 
Af H e/v showed to be satisfactory, at least in simulations per- 
formed with the Plummer softening (GBSPH). For a fixed Af 
the energy conservation deteriorates when s is reduced. If Af is 
increased and s is held fixed, an improvement in energy con- 
servation is observed. These parameters have been chosen to 
obtain an energy conservation of better than 0. 1 per cent for 
acceptable models, see Tabled 

In agreement with lTheisI (1998), the Plummer softening is 
less collisional and has a larger relaxation time than the spline 
softening, even when the correspondence is made by equating 
their potentials depths. He shows, however, that the same re- 
laxation time can be obtained with both the spline softening 
and the Plummer softening for h k. 3e. As a consequence, in 
order to maintain the same degree of relaxation and to resolve 
the vertical structure of the disc when the spline softening is 
used, many more particles are required than for the Plummer 
model. On the other hand, the Plummer softening, in compar- 
ison with the s pline softening, has a stronger stabilizing effect 
( Ro meol 1 99711 1 9981) . This may lead, as we have seen, to a late 
bar emergence. Further investigations of galaxy model stability 
using adaptive softening kernels would be worthwhile ( Dehnen 

Eoolt . 

This work was in part motivated to explain the difference 
in the characteristics of the resulting bars in various authors 
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galaxy collis ion simulations. For example, works of Barnes 
( 1992l ll998l) show very strong and prominent bars that formed 
after the first encounter. On the other hand, i n calculations 
made by Dubins ki. Mihos & Hernauisil |T996) the apparent 
bars are not observed. Although the difference may be due to 
different galaxy models, we decided to verify the numerical as- 
pects of the problem. For example, both galaxy models have 
roughly the sam e number of particles and the time-step, but e 
used by B arnes! i 19981) is six times larger than Dubinski's. 

Our collision simulations of tidally triggered bar formation 
showed that the bar properties depend on the selected gravita- 
tional softening parameter. Within the studied ranges of s and 
N, tidal bars are shorter for smaller s. Indeed, for sufficiently 
small s and N, tidal distortions can be suppressed. For exam- 
ple, we performed some additional collision simulations with 
s = 0.0005 and At = 1/1024 for N = 40 960 and N = 163 840, 
where we were able to completely suppress the bar forma- 
tion, and obtained small dense discs. Such models are however 
highly collisional and thus cannot be considered as correct. On 
the other hand, the bar's pattern speed and the bar's length vary 
weakly with s and tend to converge to a single value already for 
N ^ 491 520. Certainly, simulations with several million par- 
ticles will appro ach further convergence an d will show many 
more bar details JO'Neill & Dubinskil2003l . but an adequate s 
and At still have to be chosen. 

Confirming the results of iBerentzen et al] d2004l) . we ob- 
serve a correlation between the bar's length and the pattern 
speed. We found also, that in simulations with large e = A(Rd), 
the bars roughly conserve their length and pattern speed at least 
within 3 Gyrs, which is time between the first and the second 
collisions when a clear bar is observed. However, for smaller e 
the bar length increases and the rotational velocity reduces. 

The collision simulations have shown that a colder disc 
(large e) forms an open spiral pattern, corresponding to a late 
type galaxy, whereas the spiral arms in a hot disc (small e) 
are tightly wound. This is explained by the dispersion rela- 
tion for WKB long-branch waves if the disc's perturbation is 
weak (Toomr ell98ltlBTnnev & Tremainell987h . Another pos- 
sible explanation of tidal response differences may lie in the 
increased central density of relaxed models for small e, which 
may create the inner Lindblad resonance that prevents the feed - 
back of swing amplification mechanism (Norman et al. 1996). 

The wiggle effect observed in our simulations could imply 
that some real barred galaxies do not necessarily have trailing 
arms. A detailed study of this effect, including gas dynamics, 
would be worth. 



6. Conclusions 

In this work we have investigated the influence of the triad of 
numerical parameters (N, s, At) on the properties of an equilib- 
rium isolated galaxy model and the dynamics of two interacting 
galaxies. The main results for the galaxy equilibrium models 
can be summarized as follows: 

1 . The transfer of the disc angular momentum to the bulge and 
the halo depends on e, N, and type of softening. 



2. The convergence study indicates that for N - 491 520 and 
s = [0.001 - 0.01], the range of final values of Q changes 
only by 14 per cent, while for smaller N, the s and Q range 
broadens. Further decrease of numerical heating is ob- 
served for larger N simulations, namely, for N - 1 3 10 720, 
Q only varies 5 per cent. 

3. The final value of the Toomre stability parameter decreases 
linearly with the increasing of the Plummer softening. 

4. For the same values of N, At and e, the spline softened 
simulations produce bars earlier than those performed with 
Plummer softening. 

For the galaxy collision models it may be concluded that: 

1 . The susceptibility of galaxy models to tidal bar formation 
and its dissolution depend strongly on the chosen numerical 
parameters. Therefore, carefully chosen parameters have to 
be taken for specific simulations. 

2. For sufficiently small N, the size of the bar formed during 
an encounter, may become too long or even not appear de- 
pending on the value of s. With the increase of N artificial 
effect of the softening is reduced. 

3. The short spiral arms that are formed at the ends of the bar 
show a periodic change from trailing to leading and vice- 
versa - the wiggle. 

4. The spline softening gives smaller bars than in case of the 
Plummer softening even when they have the same potential 
depth. 

5. The properties of the merger remnant are also affected by 
the softening. The remnants of cold systems maintain their 
barred thick disc shape, whereas those of hot systems pro- 
duce rather an ellipsoid. 
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